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Abstract 

We have calculated the equation of states, the viscosity and self-diffusion coefficients, and elec- 
tronic transport coefficients of beryllium in the warm dense regime for densities from 4.0 to 6.0 
g/cm^ and temperatures from 1.0 to 10.0 eV by using quantum molecular dynamics simulations. 
The principal Hugoniot is accordant with underground nuclear explosive and high power laser 
experimental results up to ~ 20 Mbar. The calculated viscosity and self-diffusion coefficients are 
compared with the one-component plasma model, using effective charges given by the average-atom 
model. The Stokes-Einstein relationship, which presents the relationship between the viscosity and 
self-diffusion coefficients, is found to hold fairly well in the strong coupling regime. The Lorenz 
number, which is the ratio between thermal and electrical conductivities, is computed via Kubo- 
Greenwood formula and compared to the well-known Wiedemann-Pranz law in the warm dense 
region. 
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I. INTRODUCTION 



The nature of compressed matter is of considerable interest for many fields of modern 
physics, including astrophysics [it, inertial confinement fusion (ICF) {2, 3|, and other related 
fields y^. Materials under a pressure greater than a few Mbar can be drived into a strongly 
coupled, partially ionized fluid states, which is defined as the so-called warm dense matter 
(WDM). Theoretical modeling and experimental detection of the high pressure behavior of 
WDM are of great challenge and are being in extensive investigations. Among various kinds 
of WDMs, warm dense beryllium (Be) is of particular current interest. The equation of states 
(EOS) and transport properties of Be are very important in ICF, due to its appearance in 
the ablator of Deuterium- Tritium (D-T) capsule. The compressibility of the capsule, laser 



absorption, and instability growt 
thermo-physical properties of Be 



:h at the fuel-ablator interface sensitively depend on the 

a. 

The EOS of Be at shock Hugoniot up to ~ 18 Mbar have been accessed by strong shock 
waves generated by underground nuclear explosives [g, l7| • Then, it is also possible to probe 
similar pressure range in the laboratory by high intensity laser js]. Despite the success of 
these techniques in detecting wide range EOS of typical matters, one should note that after 
several decades, only seven Hugoniot points are available for Be, and more experimental 
data are desired for building successful theoretical models, such as interatomic potentials or 
chemical models currently used in SESAME EOS [9]. As another parameter in determining 
the EOS, temperature, which is difficult to be measured in nuclear explosion or high power 
laser experiments, is still needed to be clarified. Apart from the EOS, the atomic diffusion 
coefficients and fluid viscosity are key ingredients to control hydrodynamic instabilities near 
interfaces 
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12|. The electronic dynamic conductivity, from which dielectric function can 
be obtained, determines a series of interactions between laser and matters Q, Q- Due to 
these key issues to be addressed, therefore, the thermo-physical properties of Be in the warm 
dense region are highly recommended to be understood in a systematic and self-consistent 
way. 



In the present work, quantum molecular dynamics (QMD) simulations [15|, ll6|], where 
electrons are fully quantum mechanically treated by finite-temperature density functional 
theory (FT-DFT), have been introduced to study warm dense Be. The EOS are extracted 
from a series of NVT assemble sampling over different densities and temperatures, then the 
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Hugoniot curve is calculated from the Rankine-Hugoniot relation. The self-diffusion coeffi- 
cient and viscosity have been computed from the trajectory by the velocity and the stress 
tensor autocorrelation function. The dynamic conductivity, from which the DC conductivity 
and electronic thermal conductance are derived, has been obtained from Kubo-Greenwood 
formula. The rest of the paper is organized as follows: In Sec. [ITl we briefly describe the 
QMD simulations and computational method in determining the atomic transport proper- 
ties and electronic dynamic conductivity; In Sec. Illlt discussions are presented for the EOS 
and transport properties; And In Sec. HVl we get our conclusions. 



II. COMPUTATIONAL METHOD 



A. Quantum Molecular Dynamics 



Our QMD simulations employed the Vienna ab initio Simulation Package (VASP) |l7l.ll8|. 
A series of volume fixed supercells including atoms, which are repeated periodically 
throughout the space, form the elements of our calculations. After Born-Oppenheimer ap- 
proximation, electrons are quantum mechanically treated through plane-wave FT-DFT. The 
interaction between electron and ion is presented by a projector augmented wave (PAW) 
pseudopotential. The exchange-correlation functional is determined within Perdew-Burke- 
Ernzerhof generalized gradient approximation. The ions move classically according to the 
forces from the electron density and the ion-ion repulsion. The system is kept in local 
thermodynamic equilibrium with the electron and ion temperatures equal (Te = Tj). The 
electronic temperature is kept through Fermi-dirac distribution of the electronic states, and 
the ion temperature is secured through Nose- Hoover thermostat fl9 |. 

We have chosen 216 atoms in the unit cell with periodic boundary condition. A range 
of densities from p = 1.84 g/cm^ to 6.0 g/cm^ and temperatures from T = 300 K to T = 
120000 K are selected to highlight the principle Hugoniot regions. The convergence of the 
thermodynamic quantities plays an important role in the accuracy of QMD simulations. In 
the present work, a plane-wave cutoff energy of 800 eV is employed in all simulations so that 
the pressure is converged within 2%. We have also checked the convergence with respect 
to a systematic enlargement of the k-point set in the representation of the Brillouin zone. 
The correction of higher-order k points on the EOS data is slight and negligible. In the 
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molecular dynamics simulations, only T point of the Brillouin zone is included, while 4 x 
4x4 Monkhorst-Pack scheme grid points are used in the electronic structure calculations. 
The dynamic simulations have lasted 10 ~ 20 ps with time steps of 0.5 ~ 1.0 fs according to 
different conditions. For each pressure and temperature, the system is equilibrated within 
1 ~ 2 ps. The EOS data are obtained by averaging over the final 5 ps molecular dynamic 
simulations. 



B. Transport properties 

The self-diffusion coefficient D can either be calculated from the trajectory by the mean- 
square displacement 

D = ^^{m)~Rm\'), (1) 

or by the velocity autocorrelation function 

1 

^=3^ mt)-Vmdt, (2) 

where Ri is the position and Vi is the velocity of the ith nucleus. Only in the long-time 
limit, these two formulas of D are formally equivalent. Sufficient lengths of the trajectories 
have been generated to secure contributions from the velocity autocorrelation function to the 
integral is zero, and the mean mean-square displacement away from the origin consistently 
fits to a straight line. The diffusion coefficient obtained from these two approaches lie within 
1 % accuracy of each other, here, we report the results from velocity autocorrelation function. 
The viscosity 

rj = lim f](t), (3) 



t—>-oo 



has been computed from the autocorrelation function of the off-diagonal component of the 



stress tensor 



20| 



m = / {Pi2iO)Puit'))dt'. (4) 
ksT Jo 

The results are averaged from the five independent off-diagonal components of the stress 
tensor Pxyy Pyz^ Pzxj {Pxx Pyy)/'^' and {Pyy P^z)/'^- 

Different from the self-diffusion coefficient, which involves single-particle correlations and 
attains significant statistical improvement from averaging over the particles, the viscosity 
depends on the entire system and thus needs very long trajectories so as to gain statistical 



4 



accuracy. To shorten the length of the trajectory, we use empirical fits [2l| to the integrals 
of the autocorrelation functions. Thus, extrapolation of the fits to t — t- oo can more effec- 
tively determine the basic dynamical properties. Both of the D and fj have been fit to the 
functional in the form of A[l — exp(— t/r)], where A and r are free parameters. Reasonable 
approximation to the viscosity can be produced from the finite time fitting procedure, which 
also serves to damp the long-time fluctuations. 

The fractional statistical error in calculating a correlation function C for molecular- 



dynamics trajectories 22] can be given by 



where r is the correlation time of the function, and Tt^aj is the length of the trajectory. In 
the present work, we generally fitted over a time interval of [0, 4r — 5r]. 

C. Dynamic Conductivity 

The key to evaluate the electrical transport properties is the kinetic coefficients. They 
are calculated using the following Kubo-Greenwood formulation: 

^"(e) = ■^^\Mv\'<IJk')?5{^k - Cfc' - e), (6) 

k,k' 

where {ipk\v\'4'k') are the velocity matrix elements, VL is the volume of the supercell, and 
are the electronic eigenvalues. The kinetic coefficients Cij in the Chester-Thellung version 



231] are given by 



A, = (-1)^+^- j dea{e){e - ^.f-'^''\-^-^l (7) 

where /z is the chemical potential and /(e) is the Fermi-Dirac distribution function. The 
electrical conductivity a is obtained as 

o = Cu, (8) 

and electronic thermal conductivity K is 

i^ = ^(/:22-|^), (9) 

where T is the temperature. Equations ([8]) and (Q are energy- dependent, then the electrical 
conductivity and electronic thermal conductivity are obtained through extrapolating to zero 



energy. In order to get converged transport coecients, ten independent snapshots, which 
are selected during one molecular djTiamic simulation at given conditions, are selected to 
calculate electrical conductivity and electronic thermal conductivity as running averages. 

III. RESULTS AND DISCUSSION 
A. The equation of states 



TABLE L Expansion coefficients Aij for tlie internal energy E (eV/atom). 





J = 


J = l 


J = 2 


i = 


0.4369 


0.5386 


0.3149 


i = 1 


-0.7316 


1.1936 


-0.1149 


i = 2 


0.3215 


-0.1241 


0.0117 




TABLE II: Expansion coefficients Bij 


for pressure P (GPa). 




j = 


J = 1 


J = 2 


i = 


-17.9930 


10.2237 


-0.7659 


i = 1 


-62.4972 


24.4101 


0.3885 


i = 2 


38.7063 


-0.2107 


-0.0218 



Wide range EOS (p from 4.0 to 6.0 g/cm^ and temperatures of 1 ~ 10 eV) have been 
calculated according to QMD simulations, and the internal energy E (eV/atom) and pressure 
P (GPa) are fitted by expansions in terms of density (g/cm^) and temperature (eV) as 
follows: 

E = J2AjP'T\ (10) 

P = Y,B.,p'T^. (11) 

The fitted coefficient for Aijand Bi j are summarized in Tables [T] and [Tll Here, the internal 
energy E at 1.84 g/civ? and T = 300 K has been taken as zero. 

Based on the fitted EOS, the principal Hugoniot curve can be derived from the Rankine- 
Hugoniot equation, which is the locus of points in {E, P, I/)-space satisfying the condition 

{Eo - E,) + ^{Vo - Vi){Po + Pi) = 0, (12) 
6 



where the subscripts and 1 denote the initial and shocked state, respectively. This relation 
follows from conservation of mass, momentum, and energy for an isolated system compressed 
by a pusher at a constant velocity. In the canonical (NVT) ensemble in which both E and 
P are temperature dependent, the locus of states which satisfies Eq. f[T^ is the so-called 
principal Hugoniot, which describes the shock adiabat between the initial and final states. 
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FIG. 1: (Color online) Hugoniot curve computed by QMD simulations (red line) are compared 
withprevious results. Underground nuclear explosive experiments by Nellis et al. [7] and Ragan 
III \&\ are labeled as open square and triangle. High power laser results by Cauble et al. [s] are 



shown as open circles. SESAME 
lines, respectively. 



and multi-phase EOS 



24l ] are shown as dashed and dotted 



The Hugoniot curve is shown in Fig. [T], where previous theoretical and experimental 
results are also provided for comparison. In the warm dense region, the nature of the 
continuous transition from condensed matter to dense plasma remains an outstanding and 
interesting issue in high-pressure physics. One way to address this issue is to measure EOS 
in the range of 0.1 to 5 TPa. The high shock pressures required to span this range have been 
reached traditionally with strong shock waves generated by underground nuclear explosives 
0, [t], then accessed by high intensity lasers js]. Good agreement is found from Fig. Fig. [1] 
between our QMD-determined EOS and those obtained experimentally. The present QMD 



EOS indicates a smooth transition from condensed matter to plasma at temperatures from 
1.0 eV to 10.0 eV, where we do not find any signs that suggest a first order plasma phase 
transition. Theoretically, SEMASE EOS [9] also shows an overall accordance with our 
results. On the contrast, the Hugoniot curve by multi-phase EOS is soften at p ~ 6 g/cm^, 
which corresponds to the imperfect join between the mixed theoretical models 24 1. 



B. Diffusion and viscosity 




FIG. 2: (Color online) Self-diffusion coefficient and viscosity as a function of time at a density of 
5.0 g/cm'^ and different temperatures of 1.0, 5.0, and 10 eV. The fits (dashed lines) are performed 
for a sample window of [0,4T-5r]. 

We have performed ab initio quantum-mechanical simulations with the FT-DFT method 
to benchmark the dynamic properties of Be in the WDM regime. An example of the QMD 
results for the self-diffusion coefficient and viscosity of 5.0 g/cm^ at temperatures of 1.0, 
5.0, and 10.0 eV are displayed with their fits in Fig. [2l The current simulations have the 
trajectory of 10~20 ps and correlation times between 100 and 200 fs. The computed error 
lies within 10% for the viscosity. Due to the fitting procedure and extrapolation to infinite 
time a total uncertainty of ~20% is estimated by experience. Since the particle average 
gives an additional advantage, the error in the self-diffusion coefficient is less than 1%. 

Idealized model, such as one-component plasma (OCP) model, concerns the interaction 
through the Coulomb potential within a neutralizing background of electrons. A large 



amount of molecular dynamics and Monte Carlo simulations based on OCP model 
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FIG. 3: (Color online) Self-diffusion coefficient (left panel) and viscosity (right panel) as a function 
of temperature at a density of 5.0 g/cm^. Only statistical error has been considered here. 

have demonstrated that physical properties like diffusion and viscosity can be represented in 
terms of coupling coefficient, which is defined by the ratio of the potential to kinetic energy: 

where Ze is the ion charge, and a = (3/47mj)^/^ is the ion-sphere radius with rii = p/M 
the number density. A memory-function has been used by Hansen et al. [30i] to analyze the 
velocity autocorrelation function to obtain the diffusion coefficient for the classical OCP: 

= 2.95r-i-34, (14) 



with ujp 



being the ion plasma frequency. 



Bastea 25|] has performed classical molecular-dynamics simulations of the OCP and fits 
his results to the form 

AT-' + + CT, (15) 



^ = 4r-2 

njMuJnO^ 



with s = 0.878, A = 0.482, B = 0.629, and C = 0.00188. As the OCP model is restricted to 
a fully ionized plasma, determination of the ionization degree for WDM permits an extension 
of the OCP formulas to cooler systems. A reasonable choice is to replace Z in Eq. (fT3ll 
with an effective charge Z. As a consequence, we have introduced average-atom (AA) model 



3l| . which solves the Hartree-Fock-Slater equation in a self-consistent field approximation 
assuming a finite temperature. In the densities and temperatures we explore, the effective 
charge Z is 2.0, which corresponds to the case of full ionization of the 2s electrons. 



QMD and OCP results for the self-diffusion coefficient, at the density of 5.0 g/cm^, are 
shown in the left panel in Fig. [31 The tendency for the self-diffusion coefficient with respect 
to temperature is similar for QMD simulations and OCP model. However, OCP model 
predicts a larger value of the diffusion coefficient compared to QMD results (~1.5 times 
bigger). Results for the viscosity, which consists contribution from interatomic potential 
and kinetic motion of particles, are plotted in the right panel of Fig. [31 The minimum of 
viscosity along temperature can be attribute to a combined effect, that is, contribution from 
interatomic potential decreases with temperature, while contribution from kinetic motion 
increases with temperature. OCP model indicates the local minimum around 2.0 eV, and 
QMD suggests the location around 3.0 eV. 
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FIG. 4: (Color online) Examination of the Stokes-Einstein relation along the coupling parameter 
in the warm dense region. Predictions by Chisolm and Wallace [sS^ are shown as blue region and 
denoted as 'CW in the figure. The flat cyan lines show the constant values of Cse for stick and 



slip boundary conditions 



3^, 



35|. 



The Stokes-Einstein relation gives a connection between the diffusion and shear viscosity: 

Dt] 



Fse[D,v] 



SB, 



(16) 



where ^5^; is a shorthand notation for the relationship between the transport coefficients 



and is a constant. Several prescriptions 



32| for determining Cse are available. Chisolm 
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and Wallace 33| have provided an empirical value of 0.18 ± 0.02 from a theory of liquids 



near melting. On the other hand, Cse-, which has been derived based on the motion of a test 



particle through a solvent, was assumed to range from l/Gvr [34] to l/An [35!] depending on 
the limits of the slip coefficient from infinity (stick) to zero (slip). Here, we have examined 
the behavior of Be over the various regimes we explore. As shown in Fig. HI we have 
plotted the Stokes-Einstein expression Fse{,D,vi) as a function of the coupling parameter 
r using the diffusion coefficients and viscosities from QMD simulations at densities of 4 
~ 6 g/cw? and temperatures from 1 ~ 10 eV. As the expected fitting error of ~ 20% 
for the viscosity, the QMD results are bounded by the classical values of Cse from below 
(slip limit) and the Chisolm- Wallace liquid metal value from above at temperatures below 
7.0 eV (corresponding to F > 4.0). The function FsE[D,ri\ at the higher temperature 
evinces a sharp increase with temperature. The near-linear rise of the diffusion coefficient 
with temperature in the whole region basically cancels the temperature dependence of the 
denominator. As a consequence, the behavior of FsE[D,ri\ is dominated by the viscosity, 
which rises sharply at high temperatures as shown in Fig. [3l As mentioned above, this 
abrupt bend in the viscosity with temperature reflects a change from potential to kinetics 
dominated regimes. 

C. Lorenz Number 

We have also computed the Lorenz number defined as 



K 



where K and a are the thermal and electrical conductivities, respectively. 7 depends on the 



screened potential and corresponds to the scattering of the electrons [23|. In a degenerate 
{9 < 1) and coupled plasma (F > 1), L or 7 is a constant (7r^/3), reaching the ideal 
Sommerfeld number, which is the value valid for metals, and the Wiedemann-Franz law is 
recovered in an elastically interacting electron system. As temperature is very high, the 
WDM enters into a nondegenerate case (^ ^ 1 and F ^ 1), the Lorenz number reaches the 
value of kinetic matter (4 or 1.5966 depending on the e-e collisions). In the intermediate 
region, no assumptions can be found for predicting the Lorenz number and one cannot deduce 
the thermal conductivity from the electrical conductivity by using Wiedemann-Franz law. 
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FIG. 5: (Color online) Calculated Lorenz number 7 as a function of degenerate parameter 9. The 
ideal Sommerfeld number is plotted as dotted line in the figure. 

In Fig. [5] we show the behavior of 7 as a function of the degeneracy parameter 6 
(9 = ksT/Ep, with Fermi energy Ep = ^^ ^^""2^°^ ^ and electron density rig). Here, we 
should stress that in QMD simulations, the electrical conductivity a and electronic thermal 
conductance K can be directly evaluated without using any assumption of Lorenz number, 
which is highly dependent on the two non-dimensional parameters 6 and F. In the present 
warm dense regime, the Lorenz number vibrates around the Sommerfeld limit at low tem- 
peratures, and as 6 increases, a departure of the Lorenz number from the ideal value can be 
observed from Fig. [51 

IV. CONCLUSION 

In the present work, clear chains have been demonstrated in investigating the thermo- 
physical properties of warm dense Be. The EOS has been calculated through ab initio 
molecular dynamic simulations, and smooth functions have been constructed to fit the QMD 
wide range EOS data, which show good agreement with the underground nuclear explosive 
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and high pulsed laser experimental results. Based on Green-Kubo relation, the self-diffusion 
coefficients and viscosity have then been determined, and as a reference, OCP model with an 
effective charge determined from the average-atom model has also been employed. A Stokes- 
Einstein relation between the viscosity and diffusion coefficient holds the general feature of 
liquids as predicted by Chisolm and Wallace in the strong coupling region (F > 4), while, 
a sharp increase has been observed as temperature arises. The Kubo-Greenwood formula 
provides an efficient way to study the electrical conductivity and electronic heat conductance 
in the warm dense regime. Through QMD simulations we have showed that the Wiedemann- 
Pranz law is satisfied for the degenerate regime. Our present results are expected to shed 
light on the hydrodynamic modeling of target implosions in ICF design. 
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